The Annals of Applied Statistics 
2010, Vol. 4, No. 1, 320-339 
DOI: 10.1214/09-AOAS288 

© Institute of Mathematical Statistics. 2010 

CAUSAL GRAPHICAL MODELS IN SYSTEMS GENETICS: 
A UNIFIED FRAMEWORK FOR JOINT INFERENCE OF 
CAUSAL NETWORK AND GENETIC ARCHITECTURE 
FOR CORRELATED PHENOTYPES 1 

By Elias Chaibub Neto, Mark P. Keller, Alan D. Attie 
and Brian S. Yandell 

University of Wisconsin-Madison 

Causal inference approaches in systems genetics exploit quanti- 
tative trait loci (QTL) genotypes to infer causal relationships among 
phenotypes. The genetic architecture of each phenotype may be com- 
plex, and poorly estimated genetic architectures may compromise the 
inference of causal relationships among phenotypes. Existing meth- 
ods assume QTLs are known or inferred without regard to the phe- 
notype network structure. In this paper we develop a QTL-driven 
phenotype network method (QTLnet) to jointly infer a causal pheno- 
type network and associated genetic architecture for sets of correlated 
phenotypes. Randomization of alleles during meiosis and the unidirec- 
tional influence of genotype on phenotype allow the inference of QTLs 
causal to phenotypes. Causal relationships among phenotypes can be 
inferred using these QTL nodes, enabling us to distinguish among 
phenotype networks that would otherwise be distribution equivalent. 
We jointly model phenotypes and QTLs using homogeneous condi- 
tional Gaussian regression models, and we derive a graphical criterion 
for distribution equivalence. We validate the QTLnet approach in a 
simulation study. Finally, we illustrate with simulated data and a real 
example how QTLnet can be used to infer both direct and indirect 
effects of QTLs and phenotypes that co-map to a genomic region. 

1. Introduction. In the past few years it has been recognized that genet- 
ics can be used to establish causal relationships among phenotypes organized 
in networks [Schadt et al. (2005), Kulp and Jagalur (2006), Li et al. (2006), 



Received December 2008; revised July 2009. 

Supported by CNPq Brazil (ECN); NIDDK Grants DK66369, DK58037 and DK06639 
(ADA, MPK, BSY, ECN); and by NIGMS Grants PA02110 and GM069430-01A2 (BSY). 

Key words and phrases. Causal graphical models, QTL mapping, joint inference of phe- 
notype network and genetic architecture, systems genetics, homogeneous conditional Gaus- 
sian regression models, Markov chain Monte Carlo. 

This is an electronic reprint of the original article published by the 
Institute of Mathematical Statistics in The Annals of Applied Statistics. 
2010, Vol. 4, No. 1, 320-339. This reprint differs from the original in pagination 
and typographic detail. 



1 



2 



CHAIBUB NETO, KELLER, ATTIE AND YANDELL 



Chen, Emmert-Streib and Storey (2007), Liu, de la Fuente and Hoeschele 
(2008), Aten et al. (2008), Chaibub Neto et al. (2008)]. These approaches aim 
to generate a hypothesis about causal relationships among phenotypes in- 
volved in biological pathways underlying complex diseases such as diabetes. 
A key element in these methods is the identification of quantitative trait loci 
(QTLs) that are causal for each phenotype. The genetic architecture of each 
phenotype, which consists of the locations and effects of detectable QTLs, 
may be complex. Poorly estimated genetic architectures may compromise 
the inference of causal relationships among phenotypes. Existing methods 
that estimate QTLs from genome scans that ignore causal phenotypes bias 
the genetic architecture by incorrectly inferring QTLs that have indirect 
effects. 

In this paper we propose a novel framework for the joint inference of phe- 
notype network structure and genetic architecture (QTLnet). We model phe- 
notypes and QTL genotypes jointly using homogeneous conditional Gaussian 
regression (HCGR) models [Lauritzen (1996)]. The genetic architecture for 
each phenotype is inferred conditional on the phenotype network. Because 
the phenotype network structure is itself unknown, the algorithm iterates 
between updating the network structure and genetic architecture using a 
Markov chain Monte Carlo (MCMC) approach. The posterior sample of 
network structures is summarized by Bayesian model averaging. To the best 
of our knowledge, no other proposed method explicitly uses an inferred net- 
work structure among phenotypes when performing QTL mapping. Tailoring 
QTL mapping to network structure avoids the false detection of QTLs with 
indirect effects and improves phenotype network structure inference. 

We employ a causal inference framework with components of both ran- 
domized experiments and conditional probability. Randomization of alleles 
during meiosis and the unidirectional influence of genotype on phenotype al- 
low the inference of causal QTLs for phenotypes. Causal relationships among 
phenotypes can be inferred using these QTL nodes, enabling us to distin- 
guish between networks that would otherwise be distribution equivalent. 

We are particularly interested in inferring causal networks relating sets 
of phenotypes mapping to coincident genomic regions. It is widely asserted 
that alleged "hot spots" may have a "master regulator" and that most co- 
mapping is due to indirect effects [Breitling et al. (2008)]. That is, such a 
hot spot QTL could influence a single phenotype that is upstream of many 
others in a causal network; ignoring the phenotype network would result in 
a perceived hot spot. One objective of our QTLnet method is to sort out 
the direct and indirect effects of QTLs and phenotypes in such situations. 

In the next sections we develop in detail a framework for the joint infer- 
ence of causal network and genetic architecture of correlated phenotypes. 
The core idea is to learn the structure of mixed Bayesian networks com- 
posed of phenotypes (continuous variables) and QTLs (discrete variables). 
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We model the conditional distribution of phenotypes given the QTLs with 
homogeneous conditional Gaussian regression models described in Section 
2. This allows us to justify formal inference of causal direction along the 
Bayesian networks in Section 3. That is, we can reduce the size of equiv- 
alence classes of Bayesian networks of phenotypes by using driving QTL. 
In Section 4 we show how a conditional LOD score can formally measure 
conditional dependence among phenotypes and QTLs and how QTL map- 
ping can be embedded in a graphical models framework. Section 5 presents 
the MCMC approach for QTLnet. Simulation studies in Section 6 validate 
the QTLnet approach and provide an explanation for some QTL hot spots. 
Section 7 uses real data to illustrate how QTLnet can be used to infer direct 
and indirect effects of QTLs and phenotypes that co-map to a genomic re- 
gion. The Discussion puts this work in the context of open questions. Proofs 
of formal results are given in the Supplement [Chaibub Neto et al. (2009)]. 

2. HCGR genetic model. In this section we recast the genetical model 
for QTL studies as a homogeneous conditional Gaussian regression model 
that jointly models phenotypes and QTL genotypes. Conditional on the 
QTL genotypes and covariates, the phenotypes are distributed according to a 
multivariate normal distribution. The QTLs and covariates enter the HCGR 
model through the mean in a similar fashion to the seemingly unrelated 
regression model [Banerjee et al. (2008)]. However, the correlation structure 
among phenotypes is explicitly modeled according to the directed graph 
representation of the phenotype network. We derive the genetic model from 
a system of linear regression equations and show that it corresponds to a 
homogenous conditional Gaussian regression model. 

In QTL studies, the observed data consist of phenotypic trait values, y, 
and marker genotypes, m, on n individuals derived from an inbred line 
cross. Following Sen and Churchill (2001), we condition on unobserved QTL 
genotypes, q, to partition our model into genetic and recombination com- 
ponents, respectively relating phenotypes to QTLs and QTLs to observed 
markers across the genome, 

p(y, q | m) = p(y | q, m)p(q | m) = p(y | q)p(q | m) , 

where the second equality follows from conditional independence, y _IL m | q. 
That is, given the QTL genotypes, the marker genotypes provide no addi- 
tional information about the phenotypes. Estimation of the recombination 
model, p(q | m), is a well-solved problem [Broman et al. (2003)] and is not 
addressed in this paper. 

Let i = 1, . . . , n and t = 1, . . . ,T index individuals and phenotype traits, 
respectively. Let y = (yi, . . . ,y n )' represent all phenotypic trait values, y, = 
(yii, . . . , UTi)' represent the measurements of the T phenotype traits for in- 
dividual i, and let £j = (eu, . . . , £Ti)' represent the associated independent 
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normal error terms. We assume that individual i and trait t have the fol- 
lowing phenotype model: 

(2.1) y u = fJl * i + ^ PtvVvi + eu, e u ~N(0,of), 

v€pa.(yt) 

where = jjLt + Xu^t, where \i± is the overall mean for trait t, Ot is a column 
vector of all genetic effects constituting the genetic architecture of trait t 
plus any additional additive or interactive covariates, and represents 
the row vector of genetic effects predictors derived from the QTL genotypes 
along with any covariates. The notation pa(yt) represents the set of parent 
phenotype nodes of yt, that is, the set of phenotype nodes that directly affect 
Ut- Genetic effects may follow Cockerham's genetic model, but need not be 
restricted to this form [Zeng, Wang and Zou (2005)]. 

The Jacobian transformation from £j — > yj allows us to represent the 
joint density of the phenotype traits conditional on the respective genetic 
architectures as multivariate normal with the following mean vector and 
covariance matrix. 



Result 1. The conditional joint distribution of the phenotype traits 
organized according to the set of structural equations defined in (2.1) is 
y< I /<,/3,o- 2 ~ JV^n-^n -1 ), where $ = (^, . . . , fa)', (3 = {(3 tv :ve 
pa(yt),t = 1, . . . , T}, <r 2 = (erf, . . . , Oj)' , fi is the concentration matrix with 



entries given by 



1 , X-fl 



fivt fitv fisvfist 



7j is a vector with entries ^ — ^2 s ^ t ^a sl l{f _>.,,}, an d ^{t-+s} is the indicator 
function that trait t affects trait s. 

Remarks. (1) The model allows different genetic architectures for each 
phenotype. (2) The covariance structure depends exclusively in the rela- 
tionships among phenotypes since ft depends only on the partial regression 
coefficients relating phenotypes (/3's) and variances of error terms (a 2, s), 
and not on the genetic architectures defined by the 0's. (3) When the cor- 
relation between two phenotypes arises exclusively because of a pleoitropic 
QTL, conditioning on the QTL genotypes makes the phenotypes indepen- 
dent; thus, the concentration matrix of the conditional model does not de- 
pend on the genetic architecture. (4) This model can represent acyclic and 
cyclic networks. However, we focus on acyclic networks in this paper. 
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We now show that our model corresponds to a homogeneous conditional 
Gaussian regression model. The conditional Gaussian (CG) parametric fam- 
ily models the covariation of discrete and continuous random variables. Con- 
tinuous random variables conditional on discrete variables are multivariate 
normal [Lauritzen (1996)]. The joint distribution of the vectors of discrete 
(qi) and continuous (yj) variables have a density / such that 

(2.2) log/(qi,yO =g(cn) + h'(q J )y i -^K( qi ) yi /2, 

where g{<\i) is a scalar, h(qj) is a vector and K(q^) is a positive definite 
matrix. The density / depends on observed markers rtij as 

(2.3) log/(qi,yi) =logp(yi,qj | m;), 

where g(qi) = logp(qj | m,) - ±(Tlog2vr - logdet(fi) + Y%=i l4i l a t )■ Ob- 
serve that the linear coefficients h(qj) = 7^ depend on q^ through ~K t i, while 
the concentration matrix K(q,) = fi does not. Thus, our model is a ho- 
mogeneous CG model [Lauritzen (1996), page 160]. Furthermore, since our 
genetic model was derived from a set of regression equations with normal 
errors, our model is in the homogeneous conditional Gaussian regression 
parametric family. 

3. A causal framework for systems genetics. This section formalizes our 
causal inference framework for systems genetics that combines ideas from 
randomized experiments with a purely probabilistic approach to causal in- 
ference. We argue that while causal claims about the relationship between 
QTLs and phenotypes are justified by randomization of alleles during meio- 
sis and the unidirectional influence of genotype on phenotype, causal claims 
about the relationships between phenotypes follow from conditional proba- 
bility. In a nutshell, by adding QTL nodes to phenotype networks, we can 
distinguish between phenotype networks that would, otherwise, be distribu- 
tion equivalent. 

In order to formalize our approach, we first show that adding causal QTL 
nodes can break Markov-equivalence among phenotype networks by creating 
new conditional independence relationships among nodes. Second, we note 
that two models in the HCGR parametric family are distribution equiva- 
lent if and only if they are Markov equivalent. The last result together with 
Theorem 1 (see below) provide a simple graphical criterion to determine 
whether two DAGs belonging to the HCGR parametric family are distribu- 
tion equivalent. 

The analogy between the randomization of alleles during meiosis and a 
randomized experimental design was first pointed out by Li et al. (2006). 
Causality can be unambiguously inferred from a randomized experiment for 
two reasons [Dawid (2007)]: (1) the treatment to an experimental unit (geno- 
type) precedes measured outcomes (phenotypes); and (2) random allocation 
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of treatments to experimental units guarantees that other common causes 
are averaged out. The central dogma of molecular biology [Crick (1958)] sug- 
gests that QTL genotypes generally precede phenotypes. For the types of 
Eukaryotic data which we analyze (including clinical traits gene expression 
and metabolites), causality must go from QTL to phenotype. Furthermore, 
random allocation of QTL genotypes eliminates confounding from other ge- 
netic and environmental effects. 

Causal relationships among phenotypes require the additional assumption 
of conditional independence. Suppose a QTL, Q, and two phenotypes map- 
ping to Q, Y\ and Y2, have true causal relationship Q — > Y\ — > Y2. That is, 
Y2 is independent of Q given Y\. The randomization of genotypes in Q leads 
to a randomization of Y±, thus averaging out confounding affects. However, 
precedence of the 'randomized' Y\ before the 'outcome' Y2 cannot in general 
be determined a priori. Conditional independence is the key to determine 
causal order among phenotypes [Schadt et al. (2005), Li et al. (2006), Chen, 
Emmert-Streib and Storey (2007), Liu, de la Fuente and Hoeschele (2008), 
Chaibub Neto et al. (2008)]. 

In the remainder of this section we present some graphical model defini- 
tions and results needed in the formalization of our causal graphical models 
framework. A path is any unbroken, nonintersecting sequence of edges in a 
graph, which may go along or against the direction of the arrows. 

Definition 1 (d-separation). A path p is said to be d-separated (or 
blocked) by a set of nodes Z if and only if 

1 . p contains a chain i — > m — > j or a fork i <— m — > j such that the middle 
node m is in Z, or 

2. p contains an inverted fork (or collider) i—>m<—j such that the middle 
node m is not in Z and such that no descendant of m is in Z. 

A set Z is said to d-separate X from Y if and only if Z blocks every path 
from a node in X to a node in Y. X and Y are d-connected if they are not 
d-separated [Pearl (1988, 2000)]. 

Two graphs are Markov equivalent (or faithful indistinguishable) if they 
have the same set of d-separation relations [Spirtes, Glymour and Schemes 
(2000)]. The skeleton of a causal graph is the undirected graph obtained by 
replacing its arrows by undirected edges. A v-structure is composed by two 
converging arrows whose tails are not connected by an arrow. 

Theorem 1 (Detecting Markov equivalence) . Two directed acyclic graphs 
(DAGs) are Markov equivalent if and only if they have the same skeletons 
and the same set of v- structures [Verma and Pearl (1990)]. 
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Two models are likelihood equivalent if f(y | Mi) = f(y \ M 2 ) for any data 
set y, where f(y \ M) represent the prior predictive density of the data, y, 
conditional on model M [Heckerman, Geiger and Chickering (1995)]. In this 
paper we extend the definition of likelihood equivalence to predictive densi- 
ties obtained by plugging in maximum likelihood estimates in the respective 
sampling models. A closely related concept, distribution equivalence, states 
that two models are distribution equivalent if one is a reparametrization 
of the other. While likelihood equivalence is defined in terms of predictive 
densities (prior predictive density or sampling model evaluated on the max- 
imum likelihood estimates), distribution equivalence is defined in terms of 
the sampling model directly. Because of the invariance property of maximum 
likelihood estimates, distribution and likelihood equivalence are equivalent 
concepts in the frequentist setting. This is also true in the Bayesian setting 
with proper priors invariant to model reparameterizations. 

Suppose that for each pair of connected phenotypes in a graph there exists 
at least one QTL affecting one but not the other phenotype. Denote this new 
graph with QTLs included by the "extended graph." The next result shows 
that, in this particular situation, we can distinguish between causal models 
belonging to a Markov equivalent class of phenotype networks. 

Result 2. Consider a class of Markov equivalent DAGs Q. Let Y\ and 
Y 2 be any two adjacent nodes in the graphs in Q. Assume that for each such 
pair there exists at least one variable, Q, directly affecting Y\ but not Y 2 . 
Let Qe represent the class of extended graphs. Then the graphs in Qe are 
not Markov equivalent. 

As an illustrative example consider the following class of Markov equiva- 
lent models: Q = {Y x Y 2 ->• Y3, Y\ <- Y 2 <- Y 3 , Y\ <- Y 2 -> Y 3 }. These causal 
models are Markov equivalent because they have the same set of conditional 
independence relations, namely, Y\ JL Y3 \ Y 2 . In accordance with Theorem 
1, the three models have the same skeleton, Y\ — Y% — Y 3 , and the same set of 
v-structures (no v-structures). Now consider one QTL, Q, affecting Y2 but 
not Y\ and Y3. Then Qe is composed by 



Observe that these models still have the same skeleton but different sets 
of v-structures: Y\ — > Y% 4— Q, Q — > Y2 <— Y3 and 0, respectively. The next 
result guarantees that for the HCGR parametric family Markov equivalence 
implies distribution equivalence and vice-versa. 



Q 

\ 

Y 1 — Y 2 — Y 3 



Q 
\ 

Yl ^y 2 — y 3 



Q 

\ 

Yi — Y 2 — Y 3 . 
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Result 3. For the HCGR parametric family, two DAGs are distribution 
equivalent if and only if they are Markov equivalent. 

It follows from Results 2 and 3 that by extending the phenotype network 
to include QTLs we are able to reduce the size of and equivalence class of 
graphs (possibly to a single network). Furthermore, if we consider Theorem 1 
and Result 3 together, we have the following: 

Result 4. For the HCGR parametric family, two DAGs are distribution 
equivalent if and only if they have the same skeletons and same sets of v- 
structures. 

Result 4 provides a simple graphical criterion to determine whether two 
HCGR models are distribution equivalent. This allows us to determine dis- 
tribution equivalence by inspection of graph structures without the need to 
go through algebraic manipulations of joint probability distributions as in 
Chaibub Neto et al. (2008). 

4. QTL mapping and phenotype network structure. In this section we 
show that the conditional LOD score can be used as a formal measure of 
conditional independence relationships between phenotypes and QTLs. Even 
though in this paper we restrict our attention to HCGR models, conditional 
LOD profiling is a general framework for the detection of conditional inde- 
pendencies between continuous and discrete random variables and does not 
depend on the particular parametric family adopted in the modeling. Con- 
trary to partial correlations, the conditional LOD score does not require the 
assumption of multinormality of the data in order to formally test for inde- 
pendence [recall that only in the Gaussian case, a zero (partial) correlation 
implies statistical (conditional) independence], and it can handle interactive 
covariates. 

The conditional LOD score is defined as 

LOD(?/, q | x) = LOD(y, q, x) - LOD(y, x) 

(4.1) 

= log 10 

where /(•) represents a predictive density (a maximized likelihood or the 
prior predictive density in a Bayesian setting). It follows directly from this 
definition that 

(4.2) LOD(y, (? | 2; )=0 <* f(y\q,x) = f(y\x) Y JL Q | X. 

Therefore, we can use conditional LOD scores as a formal measure of inde- 
pendence between continuous (Y) and discrete (Q) random variables, con- 
ditional on any set of variables X, that could be either continuous, discrete 
or both. 



\f{y\q,x)\ (f(y\x)\ 



CAUSAL GRAPHICAL MODELS IN SYSTEMS GENETICS 



9 



Furthermore, the conditional LOD score can be used to formally test for 
conditional independence in the presence of interacting covariates (denoted 
by X ■ Q) since 

(4.3) LOD(y,q\x,xq) = log 10 { ' ' *> } - lo gl0 { I<jM } = 

if and only if Y JL {Q, X ■ Q} \ X. This is a very desirable property since, in 
general, testing for conditional independence in the presence of interactions 
is not straightforward. For example, Andrei and Kendziorski (2009) point 
that in the presence of interactions, there is no one-to-one correspondence 
between zero partial correlations and conditional independencies, even when 
we assume normality of the full conditional distributions. 

Traditional QTL mapping focuses on single trait analysis, where the net- 
work structure among the phenotypes is not taken into consideration in the 
analysis. Thus, single-trait analysis may detect QTLs that directly affect 
the phenotype under investigation, as well as QTLs with indirect effects, 
affecting phenotypes upstream to the phenotype under study. Consider, for 
example, the causal graph in Figure 1. The outputs of single trait analysis 
when Figure 1 represents the true network are given in Figure 2. 

Now let's consider QTL mapping according to the phenotype network 
structure. When the phenotype structure corresponds to the true causal 
model, we avoid detecting indirect QTLs by simply performing mapping 
analysis of the phenotypes conditional on their parents. For example, in 
Figure 1, if we perform a mapping analysis of I5 conditional on I2, I3 and 
Y4, we do not detect Q±, Q2 and Q4 because I5 _LL Q\ | 12,13,14, I5 JL 
Q 2 I 12,13,14 and Y 5 JL Q 4 | Y 2 ,Y 3 ,Y±. We only detect Q 5 since Y 5 -\IQ 5 \ 

l2,*3,l4. 

Q2 -^Y 2 

/ \ 

Fig. 1. Example network with five phenotypes and four QTLs. 
Qi Qi Qi Qi Qi Qa Qi . Qi Qa Q5 

V A > Y A > "*'• A M " 

Yi Y 2 Y s Y 4 Y 5 



Fig. 2. Output of a single trait QTL mapping analysis for the phenotypes in Figure 1. 
Dashed and pointed arrows represent direct and indirect QTL /phenotype causal relation- 
ships, respectively. 
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Now consider Figure 3(a). If we perform a mapping analysis of Y$ condi- 
tional on Y3 and Y4, we still detect Q± and Q2 as QTLs for Y5, since failing 
to condition on Y<i leaves the paths Q\ — > Y\ — > Y2 — > Yj> and Q2 — > Y2 — > I5 
in Figure 1 open. In other words, Q\ and Q2 are d-connected to I5 con- 
ditional on (13,14) in the true causal graph. Mistakenly inferring that a 
QTL has a direct effect when in reality it indirectly affects the phenotype is 
problematic, but not a serious concern. 

On the other hand, if we map an upstream phenotype conditional on 
downstream phenotypes, we could infer that downstream QTLs are causal. 
This would be a serious problem, as it would dramatically reverse the causal 
flow. Consider, for example, Figure 3(b). If we perform a mapping analysis of 
Y4 conditional on Y\ , Y3 and Y5 , we incorrectly detect Q5 as a QTL for Y4 be- 
cause in the true network the paths I4 — > Y5 <— Q5 and I4 ■<— 13 — >• I5 <— Q5 
in Figure 1 are open when we condition on Y5. That is, if we perform mapping 
analysis of a phenotype conditional on phenotypes located downstream in 
the true network, we induce dependencies between the upstream phenotype 
and QTLs affecting downstream phenotypes, and we erroneously conclude 
that a downstream QTL affects an upstream phenotype. However, a model 
with reversed causal relationships among phenotypes and incorrectly hav- 
ing downstream QTLs detected as direct QTLs for the upstream node will 
generally have a lower marginal likelihood score than the model with the 
correct causal order for the phenotypes and correct genetic architecture. 
Therefore, in practice, our model selection procedure protects against this 
type of mistake. 

5. QTLnet algorithm. In this section we propose a statistical framework 
(QTLnet) for the joint inference of phenotype network structure and genetic 
architecture in systems genetics studies. Work to date in genetical network 
reconstruction has treated the problems of QTL inference and phenotype 
network reconstruction separately, generally performing genetic architecture 
inference first, and then using QTLs to help in the determination of the 
phenotype network structure [Chaibub Neto et al. (2008), Zhu et al. (2008)]. 

(a) Q2t*-Y 3 (b) Q 2 -^Y 2 



Q x - v Y x ->- Y 3 y 5 ■<- Qs Qi - v Yi Y 3 Y 5 Q 5 

\ \ / \ I ^ 

Q4^Y 4 Qa-^Y a 

Fig. 3. QTL mapping tailored to the network structure, (a) and (b) display the re- 
sults of QTL mapping according to slightly altered network structures from Figure 1. 
Dashed, pointed and wiggled arrows represent, respectively, direct, indirect and incorrect 
QTL /phenotype causal relationships. 
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As indicated in the previous section, such strategy can incorporate QTLs 
with indirect effects into the genetic architecture of phenotypes. 

The great challenge in the reconstruction of networks is that the graph 
space grows super-exponentially with the number of nodes, so that exhaus- 
tive searches are impractical even for small networks, and heuristic ap- 
proaches are needed to efficiently traverse the graph space. The Metropolis- 
Hastings algorithm below integrates the sampling of network structures 
[Madigan and York (1995), Husmeier (2003)] and QTL mapping. 

Let M. represent the structure of a phenotype network composed of T 
nodes. The posterior probability of a specified structure is given by 

(5.1) p(.M|y,q)- 



where the marginal likelihood 

(5.2) p(y | d,M) = J r P(y I q,r,7W)p(r I M)dT 

is obtained by integrating the product of the prior and likelihood of the 
HCGR model with respect to all parameters T in the model. Assuming that 
the phenotype network is a DAG, the likelihood function factors according 
to M. as 

(5-3) p(yi | cn,T,M) = Y[p(y ti | qti,pa(y t )), 

t 

where 

(5.4) p(y ti | q ti ,pa(y t )) = AM /4 + ^ PtkVH,°t 

J/fe£pa(yt) 

and the problem factors out as a series of linear regression models. (Note 
that QTL genotypes enter the model through /z^.) 

We estimate the posterior probability in (5.1) using a Metropolis-Hastings 
algorithm detailed in Section 1 of the Supplement. The M-H proposals, 
which make single changes (add or drop an edge, or change causal direction), 
require remapping of any phenotypes that have altered sets of parent nodes. 
The accept/reject calculation involves estimation of the marginal likelihood 
conditional on the parent nodes and newly mapped QTL(s). 

Because the graph space grows rapidly with the number of phenotype 
nodes, the network structure with the highest posterior probability may 
still have a very low probability. Therefore, instead of selecting the network 
structure with the highest posterior probability, we perform Bayesian model 
averaging [Hoeting et al. (1999)] for the causal relationships between pheno- 
types and infer an averaged network. Explicitly, let A uv represent a causal 
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relationship between phenotypes u and v, that is, A uv = {Y u Y V ,Y U <— 
Y v ,Y u7 ^ Y v and Y U ^Y V }. Then 

p(A uv | y) = ^2p{A uv I M k ,y,q)p(M k | y,q) 

(5.5) * 

= 22±{A UV £ M k }p(M k \y,q). 

k 

The averaged network is constructed by putting together all causal rela- 
tionships with maximum posterior probability or with posterior probability 
above a predetermined threshold. 

6. Simulations. In this section we evaluate the performance of the QTLnet 
approach in simulation studies of a causal network with five phenotypes and 
four causal QTLs. We consider situations with strong or weak causal sig- 
nals, leading respectively to high or low phenotype correlations. We show 
that important features of the causal network can be recovered. Further, this 
simulation illustrates how an alleged hot spot could be explained by sorting 
out direct and indirect effects of the QTLs. 

We generated 1000 data sets according to Figure 1. Each simulated cross 
object [Broman et al. (2003)] had 5 phenotypes simulated for an F2 pop- 
ulation with 500 individuals. The genome had 5 chromosomes of length 
100 cM with 10 equally spaced markers per chromosome. We simulated one 
QTL per phenotype, except for phenotype Y% with no QTLs. The QTLs Qt, 
t = 1,2,4,5, were unlinked and placed at the middle marker on chromosome 
t. 

Each simulated cross object had different sampled parameter value combi- 
nations for each realization. In the strong signal simulation, we sampled the 
additive and dominance effects according to f/[0.5,l] and t7[0,0.5], respec- 
tively. The partial regression coefficients for the phenotypes were sampled 
according to j3 uv ~ 0.5C/[— 1.5, — 0.5] +0.5Z7[0.5, 1.5]. In the weak signal sim- 
ulation, we generated data sets with additive and dominance effects from 
£/[0,0.5] and £/[0, 0.25], respectively, and partial regression coefficients sam- 
pled with (3 UV ~ U[— 0.5, 0.5]. The residual phenotypic variance was fixed at 
1 in both settings. 

We first show the accuracy of the mapping analysis in our simulated data 
sets. We used interval mapping with a LOD score threshold of 5 to detect 
significant QTLs. Table 1 shows the results of both unconditional and QTL 
mapping according to the phenotype network in Figure 1. In the strong sig- 
nal setting, the unconditional mapping often detected indirect QTLs, but 
the mapping of phenotypes conditional on their parent nodes increased de- 
tection of the true genetic architectures. In the weak signal simulation, the 
unconditional mapping did not detect indirect QTLs in most cases, but we 
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Table 1 

Frequencies of QTL detection for both unconditional (top half) and conditional (bottom 

half) QTL mapping according with the true phenotype network structure in Figure 1. 
Results, for each simulation, based on 1000 simulated data sets described in the text. The 
expected architecture is the set of d-connected QTLs for the phenotype conditioning with 

respect to the network in Figure 1 



Phenotypes 




Strong signal 






Weak signal 




Expected 
architecture 


Qi 




Qa 




Qi 


Q-2 


Qa 




Yi 




0.997 


0.000 


0.000 


0.000 


0.431 


0.000 


0.000 


0.000 


{Qi} 


Yi 




0.884 


0.930 


0.000 


0.000 


0.001 


0.384 


0.000 


0.000 




Y, 




0.941 


0.000 


0.000 


0.000 


0.003 


0.000 


0.000 


0.000 


Wi} 


Ya 




0.603 


0.000 


0.690 


0.000 


0.003 


0.000 


0.370 


0.000 


{Qi,Qi} 


Yr, 




0.637 


0.321 


0.321 


0.340 


0.000 


0.000 


0.001 


0.336 


{Qi, Qi, Qa, Qb} 


Y 2 


Yi 


0.001 


0.999 


0.000 


0.000 


0.000 


0.424 


0.000 


0.000 


{Qi} 


Y 3 


\Yi 


0.000 


0.000 


0.000 


0.000 


0.000 


0.000 


0.000 


0.000 





Y 4 


Y!,Y 3 


0.000 


0.000 


0.999 


0.000 


0.000 


0.000 


0.422 


0.000 


{Qa} 


Y 5 


\Y 2 ,Y 3 ,Y i 


0.000 


0.000 


0.000 


0.999 


0.000 


0.000 


0.000 


0.415 





still observe improvement in detection of the correct genetic architecture 
when we condition on the parents. 

The expected architecture contains the d-connected QTLs when condi- 
tioning (or not) on other phenotypes as indicated in the first column of 
Table 1. For instance, Q\ and Q2 are d-connected to Y%, but only Q2 is 
d-connected to Y2 when properly conditioning on Y\. Supplementary Tables 
SI, S2, S3, S4 and S5 show the simulation results for all possible conditional 
mapping combinations. 

For each simulated data set we applied the QTLnet algorithm using simple 
interval mapping for QTL detection. The ratio of marginal likelihoods in 
the Metropolis-Hastings algorithm was computed using the BIC asymptotic 
approximation to the Bayes factor (equation 1.1 on the Supplement). We 
adopted uniform priors over network structures. We ran each Markov chain 
for 30,000 iterations, sampled a structure every 10 iterations, and discarded 
the first 300 (burnin) network structures producing posterior samples of size 
2700. Posterior probabilities for each causal relationship were computed via 
Bayesian model averaging. 

Table 2 shows the frequency, out of the 1000 simulations, the true model 
was the most probable, second most probable, etc. The results show that in 
the strong signal setting the true model got the highest posterior probability 
in most of the simulations. In the weak signal setting the range of rankings 
was very widespread. 

Table 3 shows the proportion of times that each possible causal direc- 
tion (Y u — > Y v , Y u <— Y v ) or no causal relation ({Y u -ft Y V ,Y U ft Y v }) had the 
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Table 2 

Frequencies that the posterior probability of the true model was the highest, second 
highest, etc. Results based on 1000 simulated data sets described in the text 





1st 


2nd 


3rd 


4th 


5th 


6th 


7th 


8th 


9th 


10th 


11th 


12th 


>13th 


Strong 


842 


100 


21 


11 


3 


4 


3 


2 


1 


1 


3 


1 


8 


Weak 


21 


33 


18 


19 


19 


16 


17 


15 


8 


12 


13 


5 


804 



highest posterior probability for all pairs of phenotypes. The results show 
that in the strong signal simulations, the correct causal relationships were 
recovered with high probability. The results are weaker but in the correct 
direction in the weak signal setting. 

Interestingly, single trait analysis with strong signal showed that Y§ mapped 
to Qi more frequently than to Q§ (Table 1). This result can be understood 
using a path analysis [Wright (1934)] argument. In path analysis, we decom- 
pose the correlation between two variables among all paths connecting the 
two variables in a graph. Let T> uv represent the set of all direct and indirect 
directed paths connecting u and v (a directed path is a path with all arrows 
pointing in the same direction). Then the correlation between these nodes 
can be decomposed as 

COr (?/«; Uv) = ^ ] 4 > P2u4'psP2 ' ' ' 4'vp m - 1 
P£T> UV 

(6.1) 

Table 3 

Frequencies that each possible causal relationship had the highest posterior probability 
(computed via Bayesian model averaging). Results based on 1000 simulated data sets 

described in the text 



Phenotypes 




Strong signal 






Weak signal 




— > 






— > 






(1, 2) 


0.996 


0.002 


0.002 


0.594 


0.177 


0.229 


(1,3) 


0.990 


0.002 


0.008 


0.471 


0.263 


0.266 


(1,4) 


0.990 


0.001 


0.009 


0.541 


0.196 


0.263 


(1, 5) 


0.054 


0.002 


0.944 


0.028 


0.005 


0.967 


(2, 3) 


0.016 


0.022 


0.962 


0.018 


0.017 


0.965 


(2, 4) 


0.037 


0.012 


0.951 


0.018 


0.015 


0.967 


(2, 5) 


0.997 


0.003 


0.000 


0.712 


0.075 


0.213 


(3, 4) 


0.967 


0.031 


0.002 


0.482 


0.253 


0.265 


(3, 5) 


0.997 


0.003 


0.000 


0.653 


0.116 


0.231 


(4, 5) 


0.996 


0.004 


0.000 


0.670 


0.115 


0.215 
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where 4>ij = ^ij{ vaT {yj) / var (Ui)} 1 ^ 2 1S a standardized path coefficient. As- 
suming intra-locus additivity and encoding the genotypes as 0, 1 and 2 (for 
the sake of easy computation), we have from equation (6.1) that 

cor(y 5 , qi ) = ft 1>gl (fafai + Pmfoi + Aafoi + Pufafai) j | , 



cor(y 5 ,g 5 ) =/3i,g 6 



var(g 5 ) \ 
var(y 5 ) J 



1/2 



We therefore see that if the partial regression coefficients between pheno- 
types are high, and the QTL effects j3\ Sl and f3^^ 5 and QTL variances are 
close (as in the strong signal simulation), then cor (2/5, qi) will be higher than 
cor (7/5, 55) and I5 will map to Q\ with stronger signal than to Q5. 

This result suggests a possible scenario for the appearance of eQTL hot 
spots when phenotypes are highly correlated. A set of correlated phenotypes 
may be better modeled in a causal network with one upstream phenotype 
that in turn has a causal QTL. Ignoring the phenotype network can result 
in an apparent hot spot for the correlated phenotypes. Here, all phenotypes 
detect QTL Q\ with high probability when mapped unconditionally (Table 
1). In the weak signal setting, the phenotypes map mostly to their respective 
QTLs and do not show evidence for a hot spot. No hot spot was found in ad- 
ditional simulations having strong QTL/phenotype relationships and weak 
phenotype/phenotype relations (results not shown). Thus, our QTLnet ap- 
proach can effectively explain a hot spot found with unconditional mapping 
when phenotypes show strong causal structure. 



7. Network inference for a liver hot spot. In this section we illustrate 
the application of QTLnet to a subset of gene expression data derived from a 
F2 intercross between inbred lines C3H J He J and C57BL/6J [Ghazalpour 
et al. (2006), Wang et al. (2006)]. The data set is composed of genotype 
data on 1,065 markers and liver expression data on the 3421 available tran- 
scripts from 135 female mice. Interval mapping indicates that 14 transcripts 
map to the same region on chromosome 2 with a LOD score higher than 
5.3 (permutation p- value < 0.001). Only one transcript, Pscdbp, is located 
on chromosome 2 near the hot spot locus. The 14 transcripts show a strong 
correlation structure, and the correlation structure adjusting for the peak 
marker on chromosome 2, rs3707138, is still strong (see Supplementary Ta- 
ble S6). This co- mapping suggest all transcripts are under the regulation of 
a common factor. Causal relationships among phenotypes could explain the 
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strong correlation structure that we observe, although other possibilities are 
environmental factors or a latent factor that is not included. 

We applied the QTLnet algorithm on the 129 mice that had no missing 
transcript data using Haley-Knott (1992) regression (and assuming genotyp- 
ing error rate of 0.01) for the detection of QTLs conditional on the network 
structure, and we used the BIC approximation to estimate the marginal 
likelihood ratio in the Metropolis-Hastings algorithm (equation 1.1 on the 
Supplement). We adopted uniform priors over network structures. We ran 
a Markov chain for 1,000,000 iterations and sampled a network structure 
every 100 iterations, discarding the first 1000 and basing inference on a pos- 
terior sample with 9000 network structures. Diagnostic plots and measures 
(see Section 3 on the Supplement) support the convergence of the Markov 
chain. 

We performed Bayesian model averaging and, for each of 91 possible pairs 
(Y U ,Y V ), we obtained the posterior probabilities of Y u —tY v , Y u <— Y v and 
of no direct causal connection. The results are shown in Supplementary 
Table S7. Figure 4 shows a model-averaged network. 

This network suggests a key role of 1116 in the regulation of the other 
transcripts in the network. 1116 is upstream to all other transcripts, and 




Fig. 4. Model-averaged posterior network. Arrow darkness is proportional to the poste- 
rior probability of the causal relation computed via Bayesian model averaging. For each 
pair of phenotypes, the figure displays the causal relationship (presence or absence of 
an arrow) with highest posterior probability. Light grey nodes represent QTLs and show 
their chromosome number and position in centimorgans. Riken represents the riken gene 
6530401C20ffifc. 
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is the only one directly mapping to the locus of chromosome 2. We would 
have expected the cis transcript, Pscdbp, to be the upstream phenotype in 
this network. However, the data suggests Pscdbp is causal to only two other 
transcripts and that some other genetic factor on chromosome 2 may be 
driving this pathway. This estimated QTLnet causal network provides new 
hypotheses that could be tested in future mouse experiments. 

8. Discussion. We have developed a statistical framework for causal in- 
ference in systems genetics. Causal relationships between QTLs and pheno- 
types are justified by the randomization of alleles during meiosis together 
with the unidirectional influence of genotypes on phenotypes. Causal rela- 
tionships between phenotypes follows from breakage of distribution equiva- 
lence due to QTL nodes augmenting the phenotype network. We have pro- 
posed a novel approach to jointly infer genetic architecture and causal phe- 
notype network structure using HCGR models. We argue in this paper that 
failing to properly account for phenotype network structure for mapping 
analysis can yield QTLs with indirect effects in the genetic architecture, 
which can decrease the power to detect the correct causal relationships be- 
tween phenotypes. 

Current literature in systems genetics [Chaibub Neto et al. (2008), Zhu 
et al. (2008)] has considered the problems of genetic architecture and phe- 
notype network structure inference separately. Chaibub Neto et al. (2008) 
used the PC-algorithm [Spirtes, Glymour and Scheines (2000)] to first infer 
the skeleton of the phenotype network and then use QTLs to determine the 
directions of the edges in the phenotype network. Zhu et al. (2008) recon- 
structed networks from a consensus of Bayesian phenotype networks with a 
prior distribution based on causal tests of Schadt et al. (2005). Their prior 
was computed with QTLs determined by single trait analysis. 

Liu, de la Fuente and Hoeschele (2008) presented an approach based in 
structural equation models (and applicable to species where sequence in- 
formation is available) that partially accounts for the phenotype network 
structure when selecting the QTLs to be incorporated in the network. They 
perform eQTL mapping using cis, cis- trans and irans-regulation [Doss et 
al. (2005), Kulp and Jagalur (2006)] and then use local structural models to 
identify regulator-target pairs for each eQTL. The identified relationships 
are then used to construct an encompassing directed network (EDN) with 
nodes composed by transcripts and eQTLs and arrows from (1) eQTls to 
cis-regulated target transcripts; (2) cis-regulated transcripts to cis-trans- 
regulated target transcripts; and (3) irans-regulator transcripts to target 
transcripts, and from trans-eQTL to target transcripts. The EDN defines a 
network search space for network inference with model selection based on 
penalized likelihood scores and an adaptation of Occam's window [Madigan 
and Raftery (1994)]. Their local structural models, which fit at most two 
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candidate regulators per target transcript, can include indirect eQTLs in the 
genetic architecture of target transcripts when there are multiple paths con- 
necting a cis-regulator to a cis- trans-target transcript. In other words, some 
transcripts identified as cis-regulated targets may actually be cis- trans. 

Winrow et al. (2009) rely on a (nonhomogeneous) conditional Gaussian 
model and employ the standard Metropolis-Hastings (M-H) algorithm (add, 
remove or delete edges) to search over the space of DAG structures. It differs 
from our approach in how QTL detection is coupled with the M-H proposal. 
QTLnet constructs M-H proposals for edges connecting phenotypes and 
detects QTLs conditional on the network structure. On the other hand, 
Winrow et al. first detect QTLs and then construct M-H proposals for both 
phenotypes and QTL nodes. In Section 4 of the Supplement we present a 
simulation study comparing our approach to Winrow's strategy. QTLnet 
inferred the proper genetic architecture and phenotype network structure 
at a higher rate than Winrow's approach with high signal-to-noise ratio. 
In the weak signal setting, QTLnet picks up direct causal relationship at a 
higher rate, but Winrow's is better when there is no direct causal link; on 
average, they were comparable. Although Winrow's strategy also avoids the 
detection of indirect QTLs, it is prone to miss direct QTLs when the causal 
effects among phenotypes are strong, and can perform poorly depending on 
the structure of the phenotype network (see Supplement for details). 

Current interest in the eQTL literature centers on understanding the re- 
lationships among expression traits that co-map to a genomic region. It is 
often suggested that these eQTL "hot spots" result from a master regula- 
tor affecting the expression levels of many other genes [see Breitling et al. 
(2008)]. A path analysis argument suggests that if the correlation structure 
between the phenotypes is strong because of a strong causal signal, a well- 
defined hot spot pattern will likely appear when we perform single trait 
analysis. Our simulations and real data example suggest that this is the 
situation where the QTLnet algorithm is expected to be most fruitful. 

The QTLnet approach is based on a Metropolis-Hastings algorithm that 
at each step proposes a slightly modified phenotype network and fits the 
genetic architecture conditional on this proposed network. Conditioning on 
the phenotype network structure should generally lead to a better inferred 
genetic architecture. Likewise, a better inferred genetic architecture should 
lead to a better inferred phenotype structure (models with better inferred 
genetic architectures should have higher marginal likelihood scores. A poorly 
inferred genetic architecture may compromise the marginal likelihood of a 
network with phenotype structure close to the true network). 

Because the proposal mechanism of the Metropolis-Hastings algorithm is 
based in small modification of the last accepted network (addition, deletion 
or reversion of a single edge), the mixing of the Markov chain is generally 
slow and it is necessary to run long chains and use big thinning windows 
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in order to achieve good mixing. This is a bottleneck to the scalability of 
this approach. We therefore plan to investigate more efficient versions of the 
Metropolis-Hastings algorithm for network structure inference. In particu- 
lar, a new and more extensive edge reversal move proposed by Grzegorczyk 
and Husmeier (2008) and an approach based in a Markov blanket decompo- 
sition of the network [Riggelsen (2005)]. 

Our method is currently implemented with R code. The analysis of the 
real data example took over 20 hours in a 64 bit Intel(R) Core(TM) 2 Quad 
2.66 GHz machine. We can handle up to 20 phenotype nodes, at this point. 
The run time for the code, written in R, should be substantially improved 
as we optimize code, converting key functions to C (under development). 
Nonetheless, because the number of DAGs increases super-exponentially 
with the number of phenotype nodes, scaling up the proposed approach 
to large networks will likely be a very challenging task. Current R code is 
available from the authors upon request. 

One of the most attractive features of a Bayesian framework is its ability to 
formally incorporate prior information in the analysis. Given the complexity 
of biological processes and the many limitations associated with the partial 
pictures provided by any of the "omic" data sets now available, incorporation 
of external information is highly desirable. We are currently working in the 
development of priors for network structures. 

The QTLnet approach can be seen as a method to infer causal Bayesian 
networks composed of phenotype and QTL nodes. Standard Bayesian net- 
works provide a compact representation of the conditional dependency and 
independencies associated with a joint probability distribution. The main 
criticism of a causal interpretation of such networks is that different struc- 
tures may be likelihood equivalent while representing totally different causal 
process. In other words, we can only infer a class of likelihood equivalent 
networks. We have formally shown how to break likelihood equivalence by 
incorporating causal QTLs. 

We have focussed on experimental crosses with inbred founders, as the 
recombination model and genetic architecture are relatively straightforward. 
However, this approach might be extended to outbred populations with some 
additional work. The genotypic effects are random, and the problem needs 
to be recast in terms of variance components. 

In this paper we only consider directed acyclic graphs. We point out, 
however, that the HCGR parametric family accommodates cyclic networks. 
The QTLnet approach can only infer causality among phenotypes that are 
consistent with the assumption of no latent variables and no measurement 
error. However, these complications can impact network reconstruction. We 
are currently investigating extensions of the proposed framework along these 
lines. 
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SUPPLEMENTARY MATERIAL 

Supplement to "Causal graphical models in systems genetics: A unified 
framework for joint inference of causal network and genetic architecture for 
correlated phenotypes" (DOI: 10.1214/09-AOAS288SUPP; .pdf). The Sup- 
plement article presents: (1) the Metropolis-Hastings algorithm for QTLnet; 
(2) supplementary tables for the simulations and real data example; (3) con- 
vergence diagnostics for the real data example; (4) comparison with Winrow 
et al. (2009); and (5) the proofs of Results 1, 2, 3 and 4. 
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